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We study the time evolution of two wave packets prepared at the same initial state, but evolving 
under slightly different Hamiltonians. For chaotic systems, we determine the circumstances that 
lead to an exponential decay with time of the wave packet overlap function. We show that for 
sufficiently weak perturbations, the exponential decay follows a Fermi golden rule, while by mak- 
ing the difference between the two Hamiltonians larger, the characteristic exponential decay time 
becomes the Lyapunov exponent of the classical system. We illustrate our theoretical findings by 
investigating numerically the overlap decay function of a two-dimensional dynamical system. 
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I. INTRODUCTION 

Over the last two decades the quest for quantum finger- 
prints of classical chaotic behavior has been a key subject 
of investigation in quantum chaos As a result, sig- 

natures of the classical underlying dynamics were iden- 
tified in the spectra, wave functions, and time evolution 
of a large set of quantum systems. However, one of the 
simplest indications of classical chaos, namely the Lya- 
punov exponent, remained unrelated to the quantum dy- 
namics A clear advance in this direction has been 
made recently by Jalabert and Pastawski j|, who pro- 
posed that the classical Lyapunov exponent is measured 
by the decay rate of an overlap between perturbed and 
unperturbed quantum states evolving from the same ini- 
tial state. Their work triggered several numerical studies 
[^]-^| whose results are not always in line with the original 
predictions of Ref. B . The main goal of this paper is to 
clarify the range of applicability of these predictions and 
to understand under which conditions it is possible to ex- 
tract a classical Lyapunov exponent from the quantum 
evolution of a system. 

The object of study is the comparison between the time 
evolution of a wave packet under a given system Hamil- 
tonian Ho and the corresponding evolution for a different 
Hamiltonian H = Hq + V. Formally this can be quanti- 
fied by the overlap amplitude 



0(t) = {ip\e^p{iHt/n)e^>{-iH t/fi)\ij}) 
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where, for the initial state = \ip(0)), it is chosen a the 
Gaussian wave packet 
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centered at i"o and with initial momentum po. The pur- 
pose of such parameterization is twofold: The initial 



momentum po sets the wave packet mean energy range 
at which we define (classically) the Lyapunov exponent, 
whereas the choice of a Gaussian wave packet (with finite 
width a) makes the theoretical considerations tractable 
within the semiclassical approximation. 

The amplitude overlap in Eq. (Q) can be interpreted in 
the following two different, though formally equivalent, 
ways: (a) A wave packet is prepared at the time t = 
and let to evolve under Hq till a time t > 0. The result- 
ing state is then propagated backwards in time under 
the Hamiltonian H till t = 0. Under such construction, 
|0(t)| 2 gives the return probability. This is the picture 
described by Ref. 01 which was inspired by some recent 
NMR experiments These experiments explore the 

scenario that, under certain circumstances, it is possible 
to evolve backwards in time a complex quantum system. 
This is in the spirit of the gedanken experiment at the 
origin of the Boltzmann-Loschmidt controversy [ fill and, 
for that reason, we call |0(i)| 2 the Loschmidt echo. Since 
the system is not isolated, |0(i)| 2 is expected to decay 
as t increases. The construction given by Eq. ([!]) can 
be regarded as a way to capture the physical effect of 
coupling the system to a complex time-dependent envi- 
ronment, and hence relate |0(i)| 2 to dephasing p^| , ^3[ . 
(b) Alternatively, one can regard 0(t) as the overlap am- 
plitude of an initial state \ip) propagated forward in time 
under H , with the same initial state propagated with 
H. This interpretation is closely related to the concept 
of fidelity 

Let us now state the main finding of Ref. |ij . There it 
is was shown that, after a suitable averaging (which shall 
be discussed in the foregoing section), the return proba- 
bility or fidelity can be separated into two contributions, 



M(t) = \0(t)\ 2 = M 1 (t)+M 2 (t), 



both described in the long-time limit as 



(3) 
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Mi(t) oc exp(-ait). 



(4) 



The decay rate at\ depends on the properties of the per- 
turbation V = H ~Hq, while 012 is the classical Lyapunov 
exponent associated to the dynamics of Hq, provided V 
is classically weak. Depending on V and A the decay can 
be dominated by either Mi(t) or M 2 {t). In this study we 
show under which conditions it is possible to extract A 
from the analysis of the average fidelity M(t). 

The structure of this paper is as follows. In Section 
H we describe the model we use to obtain A from the 
quantum time evolution. Section |lll| presents the analy- 
sis of the different decay processes that govern the fidelity 
M(t). There we show that Mi(t) is nothing else than the 
Fermi golden rule. The classical and quantum relevant 
scales to the problem are discussed. In particular, we 
show under which circumstances is it possible to observe 
the Lyapunov decay. The numerical results verifying a 
Lyapunov decay for our dynamical system are presented 
in Section IV . We then conclude in section ^ by relating 
our findings to the recent papers mentioned above. 



II. THE MODEL 

To investigate the dependence of the Loschmidt echo 
on the magnitude of an external perturbation we use as 
the unperturbed system the smooth "billiard" stadium 
introduced in Ref. Jl7],[l8| . This model consists of a two- 
dimensional Hamiltonian Hq — p 2 /2m + U(r) with the 
potential given by 



\2v 



U(r) = U x{ . (y/ R ) 

{[{x-dY+yiyi?) , x 



x < 0, 
< x < d, 

> d. 



(5) 



In addition, U(r) —00 whenever y < 0. The exponent 
v sets the slope of the confining potential. For v = 1 
the smooth stadium is separable and thus integrable. As 
the value of v is increased, the borders become steeper. 
In the limit of v — > 00, the stadium gains hard walls, 
becoming the well-know Bunimovich billiard, one of the 
paradigms of classical chaotic systems. (Actually, we 
consider a quarter of a stadium in order to avoid fea- 
tures related to parity symmetries. Thus, by varying 
v, we can tune the system dynamics from integrable to 
chaotic. 

In order to make the presentation more concise, we 
choose units such that U$ = 1 and m = 1/2. Thus, for 
R = d = 1 the equipotential U{x 1 y) = 1 corresponds 
to the border of the stadium with unit radius and unit 
length. For any value of the energy E the equipotential 
U(x, y) = E gives the classical turning points, defining 
the allowed area A = A(E). This area is an important 
parameter in the discussion of our numerical and analyt- 
ical results. Any exponent in the range 1 < v < 2 al- 
ready leads to a mixed phase space, i.e., a situation with 
both regular and chaotic motions present. In particular, 



for v > 2, d = 1, and total energy E = 1 the classical 
dynamics is predominantly ergodic, although small rem- 
nants of integrability still exist. These observations are 
illustrated by the Poincare surfaces of section displayed 
in Fig. | 
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FIG. 1. Poincare surface of section for the smooth stadium 
billiard for E = 1, R = d = 1, and (a) v = 1.5, (b) v = 2, and 
(c) v = 3. 

The global Lyapunov exponent A for two-dimensional 
systems can be easily computed by standard methods, 
such as that proposed by Benettin et al. |lj| . The evolu- 
tion of the classical trajectories was carried out numeri- 
cally using a symplectic algorithm (2(J . We computed the 
Lyapunov exponent for several values of v. At E = 1, A 
varies smoothly as a function of i/, as shown in Fig. I As 
expected, as v becomes very large A approaches the value 
of the Lyapunov exponent for the Bunimovich stadium 
billiard, namely Ahard = 0.86. 

The work of Ref. ||| used a Gaussian random back- 
ground potential as the perturbation that, once suddenly 
switched on, mimics the effects of external sources of ir- 
reversibility in the time evolution of a real system. Thus, 
static disorder played the role of the external perturba- 
tion V. Our strategy is essentially the same: we investi- 
gate M(t) numerically taking an ensemble average over 
different realizations of a disordered potential V(r). For 
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the later we choose a superposition of A/i independent 
Gaussian impurities, as in Refs. PU^]: 



V(r) 
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The vector Rj denotes the position of the j'th impu- 
rity. All impurities are uniformly distributed over an 
area A of the two-dimensional plane where the sta- 
dium resides, with concentration tii = Mi /A. The 
strengths Uj are Gaussian distributed and uncorrelated, 
i.e., UjUji = u 2 Sjji, with uj = 0. The impurity potential 
defined above is statistically characterized by the corre- 
lation function 



C(|r-r'|) = V(r)V{r') 



u 2 rii 
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Notice that impurity averaging yields V(v) = 0. 




FIG. 2. The Lyapunov exponent of the smooth stadium for 
E = 1 and R = d = 1 as a function of v. The circles are the 
results of our computation, while the continuous line serves 
as a guide to the eye. The dashed line corresponds to the 
billiard limit, Ah ar( j = 0.86. 



III. THEORETICAL BACKGROUND 

This section is devoted to the analysis of the time de- 
pendence of the fidelity M(t), explaining the origin of its 
different decay laws. We discuss the relation between the 
decay regimes associated to M(t) and the different time 
and perturbation strength scales of the system. These 
considerations solve the recent controversy between Lya- 
punov versus Fermi golden rule decay 

Let us start giving a more precise definition to Mi(t) 
appearing in Eq. (||), namely, 



Mx{t) = 0(t) and M 2 (t) = \0{t)\ 2 - 0{t) . (8) 

As it was already shown semiclassically Q, both Mi(t) 
and M%{t) exhibit an exponential decay in time, but dif- 
ferent in nature. We show in the sequel that the prevail- 
ing decay law is determined by the perturbation strength, 
as well as the time range under consideration. 



A. The semiclassical approximation scheme 

The best way to identify in 0(t) manifestations of the 
classical underlying dynamics is to use a semiclassical 
approximation. This is the essence of Ref. Q, which 
presents a complete calculation scheme for 0{t) in the 
case of a chaotic H and a "weak" perturbation V. The 
starting point is the Gutzwiller semiclassical propagator, 
casted in terms of a sum over all classical trajectories s 
going from r' to r in the time interval t: 



a 
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s(r,r' ;t) 



(9) 



where S denotes the action given by the integral of the 
Lagrangian, 5^(r, r';t) = J Q dt'L(V). The superscript 
V stands for the perturbation potential, its absence in- 
dicates V = 0. The Maslov index corresponding to the 
trajectory s is given by fi s and C s — \ det{d 2 S s / dr[dr j)\ 
accounts for the conservation of classical probability in 
going from the initial to the final position components 
i and j, respectively. To proceed analytically, it is nec- 
essary to restrict the calculations to a situation where 
it is possible to neglect the influence of the perturba- 
tion V in the coefficients C s . In general, the propaga- 
tor K (r, r'; t) describes the quantum evolution problem 
with great accuracy up to very long times, though shorter 
than the Heisenberg time . Since the features we are 
interested in arc manifest in a short time scale, the semi- 
classical propagator is an adequate approximation. 

We use K v (r,r';t) to propagate the wave packet 
ip(r', t — 0) given by Eq.(||) at t — up to an arbitrary 
time t. After a simple integration, one obtains 



ip V (r,t) = Vina 2 ^ Kj(r, r ; t) exp 

s(r,r ;t) 
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(10) 



where p s is defined by dS/dr'\ r / =ra = — p s . Equation 
( |T~0| ) is obtained under the assumption £ ^ a fc _1 , 
constraining the initial wave packet to be spatially con- 
centrated over a region smaller in diameter than the cor- 
relation length of the fluctuations in V(r). 
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We can now calculate the overlap Oit) as denned by 
Eq. (|ll) by writing an analytical semiclassical expression 
for {^jr (t)\ip(t)) . For times shorter than the Heisenberg 
time, this is possible through the diagonal approximation 
^ . This approximation is standard p2| and neglects con- 
tributions from pairs of trajectories which are different, 
namely, s ^ s' . The resulting expression reads 

0(t) = ^Jdv C s exp^A5. 



x exp 



s(r,r ;t) 



"(Ps -Po) 2 



where the action difference AS S is just 



AS S = - [ dt'V[q a (t% 
Jo 



(11) 



(12) 



Notice that phase difference accumulated along a trajec- 
tory s is solely due to the perturbation potential V. 

At this level, the fidelity M(t) is trivially written by 
taking the modulus squared of O(t), which implies in 
summing over pairs of trajectories s and s' taking into 
account the interference between phases, (AS*.,— AS s ')/fi- 
It is easy to check that V = leads to M{t) = 1, as ex- 
pected Q. The double sum we refer to can be split in 
two kinds of terms: (a) the diagonal ones, when the tra- 
jectories s and s' remain close to each other, and (b) the 
off-diagonal terms, corresponding to an unrelated pair of 
trajectories s and s' . In Ref. Q it was shown that after 
disorder averaging the diagonal contribution renders 



The distance in the impurity autocorrelation function C 
is r(t',t") — |q 5 (<') — <h(t")\. It is useful to change in- 
tegration variables to the center-of-mass (q + q')/2 and 
difference q — q' coordinates, with q = v^t and q' = vot' . 
For £ <C vC4, it is a good approximation to extend the 
integral over the coordinate difference to infinity. We 
can make further analytical progress if we specialize the 
discussion to hard-wall billiard systems, which are good 
approximations to our model Hamiltonian, particularly 
as v is increased. In this case, the integral over (q + q')/2 
yields L = v^t. As a result, one obtains Ef 



M{ c {t) oc exp(-a^), with 



(16) 



Notice that the Gaussian ansatz for AS S is not justified 
for very short times in the range of £/uo, which, in our 
case is of the same order as r = \CA/vq. Thus, we are 
unable to make predictions about Mi(t < r), and, conse- 
quently, about the constant factor multiplying exp(— ait) 
in Eq. (jl6|). The exponential decay can also be charac- 
terized by the typical length at which the quantum phase 
is not modified by the presence of impurities, 



U 2 Ui 



wo 



(17) 



This quantity is known as the elastic mean free path. 
Equation (Q7|) corrects a minor mistake in i given by 
Refs. namely a missing factor of 1/2. |25) In the 

sequel we show the relation between this semiclassical 
result and the stochastic theory. 



M 2 (t) oc - exp( 



-xt), 



(13) 



where A is the classical Lyapunov coefficient. In the long 
time limit, t 3> 1/A, the exponential decay dominates and 
M2(t) reduces to Eq. (Q). It is not within our scope to 
give details of this derivation, but it is worth mentioning 
that, after impurity averaging |2^,^4|], the calculations 
leading to Eq. ( [l3| ) rely solely on generic assumptions 
about the classical dynamics of H . 

The contribution to the fidelity coming from off- 
diagonal terms, Afi(i), can be computed using the im- 
purity average technique of Refs. [ 23| , p4| |. It amounts 
to computing the variance of the phase appearing in 
Eq. (pi]). Assuming that AS S are Gaussian distributed, 
which is reasonable for trajectories longer than £, one 
readily writes 



exp( T AS S ) = exp ( --^AS 2 ) , 



1 



(14) 



where, by recalling Eq. (|12|), the impurity average AS 2 
is written as 



AS 1 ? 



t pt 

dt' dt"C[r(t',t")}. 
o Jo 



(15) 



B. The random matrix approach 

The computation of Mi (t) by the statistical approach 
is a standard random-matrix result (see, for instance, 
Ref. (^(| or Appendix B of Ref. [^7j ) . A somewhat simi- 
lar calculation was also recently carried out by Mello and 
collaborators Notwithstanding, it is instructive to 

describe how this is done. The connection to the random 
matrix theory is made by the Bohigas' conjecture |||] 
and the fact that the classical dynamics of Ho is chaotic. 
Consequently, the matrix elements 



V nn > = (n\V(r)\n') 



(18) 



with respect to the eigenstates of Hq are Gaussian dis- 
tributed, regardless the form of V(r). With this in mind, 
we can calculate the averaged propagator 



K(t) 



-iHt 



/h 6(t) 



(19) 



This task is usually carried out in the energy representa- 
tion by introducing the Green function operator 



G(E) = 



1 



E + ii]- H 



with rj — > 0^ 



(20) 
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The formal expansion of G in powers of V and the rules 
for averaging over products of Gaussian distributed ma- 
trix elements give 



G — Go 



1 



1 - VGqVGo 



(21) 



where Go = (E + iT] — Hq)^ 1 . The matrix representation 
of G is particularly simple. In the eigenbasis of Hq it 
becomes 



Gnn'(E) — 



E + ir)-E n -Z n (E)' 
where E n is the n-th eigenvalue of Hq and 



(22) 



= VL-(Gv)n> = A n (E) - -T n (E), (23) 

with 



V. 



1 

nn' 



T n (E)=27rJ2V^S(E-E n ). 



(24) 



Here PV stands for principal value. The real part A^(E) 
only causes a small shift to the eigenenergy E n and will 
thus be neglected. Whenever the average matrix ele- 
ments V 2 n , show a smooth dependence on the indices 
n, it is customary to replace T n by its average value, 



T = 2irV 2 /A, 



(25) 



where A is the mean level spacing of the unperturbed 
spectrum. In most practical cases, T and A can be viewed 
as local energy averaged quantities. Hence, the average 
propagator in the time representation becomes 



K nn >{t) = S nn > exp -i 



.E n t 



Tt 
2h 



0(t). 



(26) 



It worth stressing that T arises from a nonperturbative 
scheme; nonetheless, it is usually associated to the Fermi 
golden rule due to its structure. 

The average propagator obtained in Eq. (|2^) is easily 
related to M\(t) by calculating (tp\K\ip). This step gives 
us also a more precise meaning to the smooth energy de- 
pendence of T(E): In our construction the latter has to 
change little in the energy window corresponding to the 
energy uncertainty of ^(r, t), which is determined by a. 
Thus, the RMT final expression for Mi(i) is 



M-f (f) = exp(— rt/7i), 



(27) 



with r given by Eq. (|25|). Equation (|27j) does not hold 
for very short times, since we neglected the smooth en- 
ergy variations of T n and A„. It is beyond the scope 
of RMT to remedy this situation, since for that purpose 



nonuniversal features of the model have to be accounted 
for. 

Despite sharing the same formal structure, it remains 
to be shown that both scmiclassical and random model 
theory are strictly equivalent. This is what we do next 
by deriving an expression for the Fermi golden rule in 
terms of the classical quantities used in Eq. ( |l6| ) . 

For chaotic systems, we can calculate the average off- 
diagonal perturbation matrix elements using the univer- 
sal autocorrelation function of eigcnstates first conjec- 
tured by Berry Q. For two-dimensional billiards this 
function reads 

(MnWnfa)) = ^MK\ri - r 2 |), (28) 

where Jq(x) is the Bessel function of zero order, k n is 
the wave number associated to the nth eigenstate of 
Ho, and A is the billiard area. Here (• • •) can be re- 
garded as the average ?/> n (ri )■)/>„ (r 2 ) obtained by sweep- 
ing R = (r 2 + i"i)/2 over a region containing several de 
Broglie wave lengths. Equivalently, one could also av- 
erage over a large number of levels, provided that k n 
does not change much on that interval. For a rigorous 
discussion on the validity of Eq. fl2g| ) and the different 
averaging procedures, see Ref. Q |. 

Recalling Eq. flTq), we can write the off-diagonal 
squared matrix elements averaged over the impurity re- 
alizations as 



Vnn'= l d r i M r 2 ^„(ri)V>„(r 2 )Vv(ri)V> 



xV(n)F(r 2 ) 



(29) 



By changing variables to R = (r 2 + i"i)/2 and r = r 2 — i"i 
and with the help of Berry's conjecture, it is straightfor- 
ward to reduce the integral in Eq. j29"|) to 



VL- = ^ l<t 2 rJ s Ak,,r)J„(k,,r\Cu). CUD 



The correlation function G is given by Eq. ( ]7|). For a 
sufficiently large billiard, £ <C A 1 / 2 , we obtain []32| 



V 2 , 

nn 



TliU 



A 



< k ~+wi (2k n k n ,e) , (3i) 



where Io(x) is the modified Bessel function of the first 
kind. For high-energy eigenstates, such that fc„£ ^> 1, 
and for states within an energy window corresponding to 
a (k n « k n i), the above expression is further simplified 
to 



riiU 



1 



A 2^Hkni 



(32) 



We can now insert V 2 n , into the left-hand side of Eq. 
(p5|). Recalling that the mean level spacing for a two- 
dimensional billiard is A = 2nh 2 /(Am) and using hk = 
mvo, we obtain 
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r u rii 



(33) 



its classical Lyapunov exponent A and the perturbation 
strength u. The results are presented in the next section. 



This is exactly the same decay rate of Eq. (|l^) . It also 
agrees with the quantum diagrammatic perturbation the- 
ory for the bulk in the disordered model. |24| 



C. Fermi golden rule and Lyapunov decay 

By employing the semiclassical approach we were able 
to address in detail two very distinct regimes of M(t). 
Such approximation is the most appropriate tool to study 
M{t) provided two conditions are met: V is (a) classi- 
cally weak, in the sense that classical perturbation theory 
applies, and (b) quantum mechanically strong, meaning 
that one can treat the actions in Eq. (|l4|) as Gaussian 
variables. In such cases, for t ^> A -1 , it was found that: 
(a) log M2(t) cx —At independently of the strength of the 
perturbation and (b) log Mi (t) cx —Ft, where T cx u 2 . 
As one varies the perturbation strength it, M(t) is domi- 
nated by the smallest of A and T. In other words, within 
the semiclassical regime, for small values of u the Fermi 
Golden Rule applies. The dependence of M{t) crosses 
over to logM(i) cx —At when T > A. Equations ( |l3| ) and 
( |l6| ) predict for which value of u this transition occurs. 

It remains to be discussed what happens to M(t) when 
u does not meet neither (a) nor (b), namely either u is 
below the Fermi Golden Rule regime or u is in the very 
opposite limit of strong perturbations, above the Lya- 
punov regime. 

Let us hrst discuss the limit of "extremely" weak per- 
turbations, where V neither significa ntly mixed the states 
of H , nor causes level crossings Here, M(t) can 

be obtained by standard quantum perturbation theory. 
This limit was studied a long time ago by Peres Jl3] , who 
found a Gaussian decay, namely log M{t) cx —(ut) 2 . It 
turns out, as illustrated by our numerical study, that this 
limit is very hard to observe, since for very short times 
M(t) decays as t 2 in all cases. 

At the opposite end there is the case of "strong" per- 
turbations, for which classical perturbation theory breaks 
down. As it is augmented the Lyapunov exponents of Hq 
and H become increasingly different. Lacking a theoret- 
ical understanding for this regime, we can only speculate 
that M(t) decays faster than in the Lyapunov regime. 
Here M(t) will strongly depend on the specific details of 
V(r). 

Figure |3| summarizes the principal predictions of this 
section. The main feature of this diagram is the plateau 
in — log M(t) versus u, characterizing the Lyapunov 
regime. For a given specific system we can predict where 
the plateau starts at low values of u. To use a quantum 
system to measure the Lyapunov exponent, it is crucial 
to know where it ends, and classical perturbation the- 
ory breaks down. For that purpose, numerical simula- 
tions were performed for the smooth stadium by varying 



o 




logu 

FIG. 3. Sketch of the expected behavior for log M(t) as 
a function of the perturbation strength u for a fixed value 
of t. The shaded fields indicate the regimes of u where the 
semiclassical approach fails. 



IV. NUMERICAL RESULTS 

In this section we present a numerical study of M(t) 
for the smooth stadium model with Gaussian impurity 
disorder introduced in Section |l|. Before showing the re- 
sults, however, we describe some technical details about 
the numerical method employed in the simulations. 

The quantum evolution of wave packets, as defined 
by Eq. (||), was carried out through the fourth-order 
Trotter-Suzuki algorithm |35|] . It is worth noticing that a 
more straightforward approach, based on a matrix repre- 
sentation of the evolution operator in terms of the eigen- 
vectors of Ho would be far less efficient. 

The method does not require spatial discretization of 
the system. However, the basis has to be such that the 
Hamiltonian matrix elements needs to involve only short 
term interactions. We thus found useful to work on a 
lattice and to represent the kinetic energy with a nearest- 
neighbor hopping term. Within the energy range we ex- 
plored, we found that a two-dimensional lattice of area 
2.1R x 1.1R provided very accurate results when we em- 
ployed N = 180 sites per unit distance R (with the in- 
tersite distance given by a — R/N), corresponding to a 
total number of 378 x 198 lattice sites. 

The range of parameter values explored in our sim- 
ulations was limited by computational cost. Moreover, 
our choice of parameters was guided by the const raints 
imposed by the semiclassical calculations of Sec. [II A. 



First, in order to include a large number of randomly 
located impurities, their correlation width £ had to be 
taken much smaller than R. Second, the semiclassical 
regime where Eq. (||) applies requires £ to be larger than 
the wave packet width a, which, in turn, has to be much 
larger than the particle wavelength A f ■ Other constraints 
arise from finite size effects. For instance, the large-time 
saturation value of the Loschmidt echo, M(t — ► oo), de- 



G 



pends on the ratio a/N. Thus, for a fixed N, it is neces- 
sary to make a as small as possible in order to guarantee 
a small value for M(t — > oo). In addition, let us recall 
that the energy spectrum of the (open boundaries) dis- 
cretized system is given by 

^-2 

£k = n o [cos(fe x a) + cos(fc v a)l . (34) 

ma 2 ma 2 

Therefore, we can only accurately recover the disper- 
sion relation of the free particle, E^ = h 2 k 2 /2m, when 
ka -C 1. All these constraints are summarized by the 
inequalities 

a<A F <er<£<i?. (35) 

The compromise between a good accuracy and a feasi- 
ble simulation time led us to set £ = 0.25.R, a = 0.18R, 
Xf = 0.07 R, and N — 180. This choice, combined with 
the values of the classical model parameters, m = 1/2 
and E = 1, gave raise to units such that Ti = O.Olli?. 
Thus, the inequalities of Eq. (|3^) were approximately 
observed in our simulations. For the quantum evolution, 
a time step St — 2ma 2 /10h — 2.8 x 10~ 4 E/h proved to 
be sufficiently small. 

It is important to make a few remarks about the av- 
eraging procedure. In the simulations, besides averaging 
over impurity configurations, we also found important to 
average over initial positions ro and directions po- The 
main reason is that numerical simulations of billiards deal 
with relatively small, confined systems and directionality 
has a strong influence in the short-time dynamics. 

The initial conditions for the quantum evolution were 
chosen from a subset that also minimized finite-size ef- 
fects. That is, we chose initial conditions that allowed for 
the observation of an exponential decay before the satu- 
ration time. For that purpose, we took 0.57? < xq < R, 
0.2R < yQ < 0.5R, and initial momentum po such that 
the first collision with the boundary occurred at x > R, 
avoiding trajectories close to bouncing ball-like modes 
along y. (Such trajectories were found to lead to strong 
non-exponential decays in M(t) for time intervals shorter 
than the saturation time.) 

In Fig. we show M(t) for v = 1.5, 2, and 3 for differ- 
ent values of the perturbation strength. In all graphs we 
see that the asymptotic decays are approximately expo- 
nential within a certain ranges of u, as predicted in Ref. 

In order to obtain the characteristic decay times, we 
fitted lnM(t) to the function ln[Aexp(-t/T^)/t + MJ. 
The fit was performed for times t > R/v, where v = 
v / 2E/m = 2 is the wave packet velocity, to exclude the 
initial, non-universal (and non-exponential) time evolu- 
tion. It is worth noticing that the usual nonlinear fitting 
procedures are rather insensitive to certain combinations 
of parameters and A. Thus, while the parameter 
could be fixed by averaging the long-time tail of the data, 
we avoided the uncertainty in A and by fixing the value 
of the fitted curve at the initial point to be exactly equal 



to the respective data value. We checked that such proce- 
dure yield values for A proportional to u~ 2 , as expected. 




2 4 6 8 10 12 14 



t [Rm/p o ] 

FIG. 4. M(t) for v = 1.5 (a), 2 (b), and 3 (c) for different 
values of the perturbation strength: u — 0.002, 0.005, 0.01, 
0.02, 0.03, 0.04, 0.05, and 0.06. 



The typical number of samples used in the averag- 
ing procedure (for each trace of the M[t) shown) was 
in the range 80-100. In fact, we observed that the num- 
ber of samples needed to obtain comparable statistical 
mean squares fluctuation for M(t) scaled with the per- 
turbation strength u. That is, the larger the perturba- 
tion, the larger the fluctuations in M(t) were. This fact 
set another practical limit to the range of perturbation 
strengths u we could investigate in our numerical simu- 
lations. 



In Fig. |5j we show the fidelity curves for the same 
perturbation strength, but different steepness of the con- 
fining potential. Notice that the fluctuations around the 
(exponential) fitted curve increase as the billiard walls 
become softer. 
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behaviors for 1/t$ at small and large values of u. 

The most pronounced feature shared by all data sets 
is the existence of a plateau around the classical Lya- 
punov exponent A, as expected. The semiclassical theory 
H predicts that this saturation should appear when the 
perturbation is quantum mechanically strong, but clas- 
sically weak. This condition, already presented in Eq. 
(35|), can be translated into the inequality A <C v/l. In- 
deed, the results of the simulations, as presented in Fig. 

, are consistent with the existence of a plateau in 
for u within this range. For weak perturbations, the data 
is also consistent with the quadratic behavior of ot\. 



V. CONCLUSIONS 







FIG. 5. The fidelity as a function of time, u — 0.01 for 
v = 1.5, 2, and 3. The number of samples used in the averag- 
ing behind the v = 1.5 curve was 80. 100 samples were used 
in 
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FIG. 6. The characteristic decay rates obtained from Fig. 
|as a function of perturbation strength. The solid curves 
correspond to the phenomenological expression, Eq. (fj6|). 

In Fig. H we have plotted the inverse characteristic de- 
cay times I/t^ obtained in the fittings as a function of 
the impurity strengths u for the three values of v. For 
comparison, we plotted the phenomenological curve 



Tphenom (^) 



1 

A 



a\(u) ' 



(36) 



where A is the classical Lyapunov (u independent) and 
a.\ = vq/ 1 is the characteristic decay rate obtained in 
Sec. [II A. Such curve matches the expected asymptotic 



We studied the time evolution of two wave packets 
prepared at the same initial state and time, but evolv- 
ing under slightly different Hamiltonians, namely Ho and 
H — Hq + V. For those systems for which the Hamil- 
tonian Ho is classically chaotic, the wave packet overlap 
decays exponentially in time, according to the semiclassi- 
cal theory [Q. For the model Hamiltonian introduced in 
Sec. H we numerically verified that the exponential decay 
is indeed observed for a broad range of typical strengths 
of V. 

Within the regime of perturbations which are quantum 
mechanically strong, but classically weak, the semiclassi- 
cal theory predicts two decay laws Q . While the first one 
is governed by the mean free path and the wave packet 
mean velocity, a± — Vo/i, the second decay law is char- 
acterized by the Lyapunov exponent, namely = A (22). 
By estimating the variance of V nn i we showed that ol\ is 
nothing else than the Fermi Golden Rule of Ref. . Our 
analytical findings are in quantitative agreement with the 
numerical results obtained from the smooth stadium. 

Finally, for sufficiently long times we were able to qual- 
itatively understand the behavior of the fidelity M(t) as 
a function of the strength u of V. For very weak u, quan- 
tum perturbation theory applies and log M(t) oc — (ut) 2 
|]l4j| . Increasing u, one enters in a regime where although 
quantum perturbation theory breaks down, the classical 
one still holds. We call this the semiclassical regime. 
Here, we quantified the crossover from the Fermi Golden 
Rule decay to the Lyapunov decay. Finally, by further 
increasing it, classical perturbation is no longer valid and 
the semiclassical approximation ceases to be useful. This 
picture is illustrated by Fig. ^| and nicely numerically ver- 
ified by Fig. ^. The plateaus obtained in the simulations 
show that it is possible to measure the classical Lyapunov 
exponent with quantum mechanics over a broad range of 
perturbation strengths. 
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